# Figure B.5: Accuracy for Different Values of rho


# 1. Load Packages ----

library(tidyverse)
library(latex2exp)

# 2. Read in Data ----

load(file = "df_voxit_individual_allrefs.RData")

# 3. Data Preparation ----

# 3.a Define Functions ----
SimulatedTurnout <- function(data,p=1,c,f,rho){
  # This function creates predicted turnout based on the model by Foellmi, Heim, 
  # and Schmid (2025).
  
  data_i <- data %>%
    group_by(id, anr) %>%
    mutate(u_rho=max(p*salience+c*difficulty,0)^(1/rho)) %>%
    ungroup() %>%
    group_by(id, datum_merge) %>%
    summarize(u_sum=sum(u_rho)^(1/(1/rho))+f,
              turnout=max(turnout)) %>%
    ungroup() %>%
    mutate(turnout_pred=ifelse(u_sum>0,1,0)) 
  
  return(data_i)  
}

EvaluateTurnout <- function(data){
  # This function takes the predicted turnout and real turnout data as input 
  # and then calculates the share of correct predictions. 
  
  df <- data %>% 
    mutate(correct=ifelse(turnout==turnout_pred,1,0))  %>%
    group_by(correct) %>%
    summarize(n=n())
  
  total <- sum(df$n)
  
  return(df$n[df$correct==1]/total)
}

# 3.b Loop ----
# Loop over all possible cases/combinations
Gridsearch <- function(data, p=1,
                       cspace=c(0:10),
                       fspace=c(0:10),
                       rhospace=c(0:10)){
  
  out <- data.frame()
  for (c in cspace){
    for (f in fspace){
      for (rho in rhospace){
        dof <- SimulatedTurnout(data=data,p=p,c=c,f=f,rho=rho)
        out <- bind_rows(out,data.frame(p=p,c=c,f=f,rho=rho,
                                        share_correct=EvaluateTurnout(dof))) 
      }  
    }  
  }  
  return(out)  
}



# 3.c Predictions ----
GetPredictions <- function(rho, data, c, f, p=1){
  # This function predicts turnout and computes the accuracy 
  
  df <- SimulatedTurnout(data=data, c=c, rho=rho, p=p, f=f) %>%
    mutate(category = case_when(
      turnout == 1 & turnout_pred == 1 ~ "true positive",
      turnout == 0 & turnout_pred == 0 ~ "true negatives",
      turnout == 0 & turnout_pred == 1 ~ "false positive",
      turnout == 1 & turnout_pred == 0 ~ "false negative"),
      cat_tf = case_when(
        turnout == 1 & turnout_pred == 1 | turnout == 0 & turnout_pred == 0 ~ "true",
        turnout == 1 & turnout_pred == 0 | turnout == 0 & turnout_pred == 1 ~ "false"))
  
  df_summary <- df %>%
    filter(!is.na(turnout)) %>%
    group_by(category, cat_tf) %>%
    summarise(nobs = n()) %>%
    ungroup() %>%
    mutate(share = nobs/sum(nobs)*100)
  df_summary$share_tf <- sum(df_summary$share[df_summary$cat_tf=="true"], na.rm = T)
  
  colnames(df_summary) <- c("category", "cat_tf", paste0("nobs:", rho), paste0("share:", rho),
                            paste0("share_tf:", rho))
  colnames(df) <- c("id", "datum_merge", paste0("u_sum:", rho), "turnout", 
                    paste0("turnout_pred:", rho), paste0("category:", rho))
  return(list(df, df_summary))
} 

GetPredictionsAll <- function(p=1, rho, c, f, data){
  # This function applies the function GetPredictions to a vector
  
  dfs <- lapply(X = rho, FUN = GetPredictions , 
                data = data, c = c, f=f)
  df_all <- Reduce(function(...) merge(..., by = c("id", "datum_merge", "turnout")), 
                   sapply(dfs, `[`,1))
  df_summary_all <- Reduce(function(...) merge(..., by = c("category", "cat_tf")), 
                           sapply(dfs, `[`,2))
  return(list(df_all, df_summary_all))
}

GetAccPlot <- function(p=1, rho, c, f, data){
  # This function applies the function GetPredictionsAll and reshapes the data
  
  tezt <- GetPredictionsAll(p=p, rho = rho, data = data, c = c, f = f)
  df_acc <- tezt[[2]] %>%
    pivot_longer(data =., cols = c(-category, -cat_tf), 
                 names_to = c(".value", "rho"), names_sep = ":") %>%
    mutate(parameters = paste0("c=", c))
  return(df_acc)
}


# 3.d Generate Data and Plot it  ---
rho <- c(0.002, 0.0025, 0.003, 0.004, 0.006, 0.008, 0.01, 
         0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1)
df_plot_all <- data.frame()
for(i in rho){
  df_plot <- GetAccPlot(rho = i, 
                        c = -0.86333197, p=1, f = f,
                        data = voxit_complete_salience[voxit_complete_salience$prop_voxit>1,])
  df_plot_all <- rbind(df_plot_all, df_plot)
}


# 4. Plot ----
ggplot(data = df_plot_all, 
       mapping = aes(x = as.numeric(rho), y = share_tf)) +
  geom_point() +
  geom_line() +
  scale_x_continuous(name = TeX("$\\rho$"),
                     limits = c(0, 1),
                     breaks = seq(0, 1, 0.1),
                     labels = seq(0, 1, 0.1)) +
  scale_y_continuous(name = "Accuracy",
                     breaks = seq(74.0,75.6, 0.2),
                     labels = paste0(seq(74.0,75.6, 0.2),"%"),
                     limits = c(73.8, 75.6)) +
  theme_bw(base_size = 11)

# 5. Save Plot ----
ggsave(filename = "FigureB5.pdf", 
       plot = last_plot(),
       width = 7.55, height = 5)


